DYNAMICS OF POINT JOSEPHSON JUNCTIONS IN A 
MICROSTRIP LINE 
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Abstract. We analyze a new long wave model describing the electrodynamics of an array of 
point Josephson junctions in a superconducting cavity. It consists in a wave equation with Dirac 
delta function sine nonlinearities. We introduce an adapted spectral problem whose spectrum gives 
the resonances in the current-voltage characteristic curve of any array. Using the associated inner 
product and eigenmodes, we establish that at the resonances the solution is described by two simple 
ordinary differential equations. 
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1. Introduction. The macroscopic state of a superconductor is described by a 
complex field the order parameter. For low T c superconductors it can be assumed 
that only the phase of the order parameter varies. The coupling of two such super- 
conductors across a thin oxide layer is described by the Josephson equations pQ. 

(1.1) V=-, 7 = S J cS in(-), 

where $ is the phase difference between the top and bottom superconductor, V and I 
are, respectively, the voltage and current across the barrier, s is the contact surface, 
J c is the critical current density and $o = h/2e is the reduced flux quantum. These 
two Josephson relations together with Maxwell's equations imply the modulation of 
DC current by an external magnetic field in the static regime and the conversion 
of AC current into microwave radiation [2j [3]. Such Josephson junctions are then 
unique electronic systems for applications like the detection of magnetic fields, ultra 
fast electronics 13! and microwave sources and signal mixers [1],[5]. 

For the applications the devices are often associated to form arrays. The junctions 
can be in parallel or in series. The series arrays can lead to synchronization^ . and 
deliver more output power for some applications. Their description is however more 
complex and we will not consider it here. Parallel arrays where the junctions are 
embedded between two superconducting planes are now relatively easy to prepare and 
the junction is protected from the atmosphere. In addition one can easily prepare an 
array with junctions of specific sizes and positions. Such non uniform arrays have 
been produced and analyzed in particular by Salez and co-workers at the Observatory 
of Paris. For these systems, the phase difference $ satisfies an inhomogeneous 2D 
damped driven sine Gordon equation [7] resulting from Maxwell's equations and the 
Josephson constitutive relations (|l.lj) . The damping is due to the normal electrons 
and the driving through the boundary conditions with an external current or magnetic 
field applied to the device. 

To model such arrays authors have used lumped models where the spatial de- 
pendence between the junctions is omitted. This obliterates the wave features of the 



•LABORATOIRE DE MATHEMATIQUES, INSA DE ROUEN, B.P. 8, AVENUE 
DE L'UNIVERSITE 76801 SAINT-ETIENNE DU ROUVRAY, FRANCE. E-MAIL: 
CAPUTO@INSA-ROUEN.FR, LOUKITCH@INSA-ROUEN.FR 

1 



2 



J.G. Caputo and L. Loukitch 



solution and does not describe well the experiments. Solving numerically the full two 
dimensional problem is of course possible, however it does not lead to understand 
simply the role of the parameters. Similar difficulties occur with global (hard) analy- 
sis. Consider for example the problem of finding the maximum current giving a static 
solution for a given magnetic field. Using such global analysis we obtained bounds [8] 
on the gradient of the solution that were independent of the area of the junctions so 
that little information could be obtained from them. To overcome these difficulties 
we recently introduced a continuous/discrete model that preserves the continuity of 
the phase and its normal gradient across the junction interface and where the phase is 
assumed constant in the junctions. The relative simplicity of the model allowed an un- 
precedented understanding of the static problem^ and gave excellent agreement with 
the complex static response of the array[10 . Additionally the model allows to solve 
the inverse problem of building a device that produces a given static behavior 

The dynamic behavior of Josephson junctions is characterized by the current- 
voltage (I-V) characteristic curve. To understand this one needs to analyze periodic 
solutions of the problem. For a homogeneous long junction Kulik[12 developed a 
formalism using a high voltage ansatz and obtained average equations describing the I- 
V curve. This approach was extended by Cirillo et al [13] who showed that a magnetic 
field r reinforces the cavity modes such that V = mr/l. Using this approach these 
authors obtained excellent quantitative agreement with their experimental results. 
For arrays of equidistant junctions, a recent study by Pfeiffer et al[14] analyzed the 
fine features of the first resonant step in terms of Cerenkov radiation between a sine- 
Gordon discrete kink (fiuxon) and a cavity mode. For one junction in a cavity, our 
theoretical study [TS] revealed that the junction could stop waves across the cavity 
or enhance them throughout. We also found kink like solutions for the problem Q]5] 
and explained some features of the current-voltage characteristics. Here we analyze 
theoretically and numerically the model in particular when there is a capacity miss- 
match between the junctions and the cavity. This capacity ratio is usually large in 
experiments because the oxide layer in the junction is about 10 Angstroms while it is 
about 0.2 Micron in the strip. Taking this miss-match into account we introduce an 
associated linear problem which enables us to predict the position of the resonances 
in the current-voltage curve for any array. This linear problem defines eigenvalues 
and eigenmodes orthogonal with respect to an inner product that we establish. At 
these resonances the solution just contains the Goldstone mode and the corresponding 
eigenmode so the dynamics is described by two simple amplitude equations that we 
present and analyze. 

The article is organized as such. After introducing the model in section 2, we analyze 
it in section 3, establish the periodicity of the current voltage curve as a function of the 
magnetic field and simplify the model using the time averaged (high voltage) solution. 
The resonances are studied in section 4 where we define the appropriate spectral 
problem for any given array and find its spectrum and associated inner product. 
Using the latter we analyze numerically the current voltage curves in section 5. In 
particular we establish and analyze the amplitude equations describing the system at 
resonance and discuss the situation for arrays containing many junctions. 

2. The model. 

2.1. The 2D problem. The device we model (see Fig J2.1|) is composed of two 
overlapping superconducting layers, in which stand smal0 Josephson junctions. Using 
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Fig. 2.1. The left panel shows the top view of a superconducting microstrip line containing 
three Josephson junctions. The parameters H, I and <j> are respectively the applied magnetic field, 
the current and the phase difference between the two superconducting layers. The right panel shows 
the associated 2D domain of size I X w containing n = 3 junctions placed at the positions y = w/2 
and x = a;, i = 1, n. 



the Josephson constitutive equations f eg J 1 . 1 [) and Maxwell's equation, one obtains the 
following inhomogeneous sine-Gordon equation for $ in the junction region Qj [2J [3] 



(2.1) 



J r sin | — 
$ 



-$t = 
R 



where Cj is the capacity of the junction per unit area, R the resistance per unit 
area due to normal electrons. The branch inductance L — (IqcIq involves the magnetic 
thickness do, a quantity used by experimentalists. In the microstrip Q — Qj, the 
Josephson and quasi-particle currents are absent so one obtains the wave equation, 



(2.2) 



QQu - -^A$ = 
^i 



where the I subscripts indicate that we are in the linear region (i.e.: outside the 
junction). For inhomogeneous circuits like the one of Fig. 12.11 these two equations 
can be written as |T51 [TB] 



(2.3) C&tt - V 



1 



L(x) 



g{x,y) 







(*o) 


+t] 



= 0, 



where g — 1 in the junctions and outside. This formulation guarantees the continuity 
of the phase and its normal gradient across the interfaces. At this point we will assume 
the same surface inductance L in the junctions and linear region. This simplifies 
greatly the formulation and can be realized in practical situations. To normalize the 
equation, we introduce the units of length and time, respectively the Josephson length 
Aj and plasma frequency u> p 



(2.4) 



A., 



1 



JcX'l Up V Jc 

Finally we normalize space, time and phase as 

(2.5) x = x/Xj, y = y/Xj, t = tui p , ^ = $/$ , 

to get the normalized 2D inhomogeneous perturbed sine-Gordon equation 



(2.6) 



fit - Ai p + v){k<p% + ai ft + sin v>) = 
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where the coefficients a and k are 



The boundary conditions are 

(2.8) ^| (fc0)=J ff+(i_ i ,)_L, ^| ( . = - =j ff_( 1 _ 2y )_L i 

(2.9) ¥>y 1(5=0) = </>yl(£=«) = 
where 

(2.10) H = H^ 1 i = I ^ = LL^ 

and I and u> are normalized by Xj. After this section all tildes will be omitted for 
simplicity. 

2.2. The ID model. This equation is difficult to analyze and its solutions can 
only be obtained numerically. In addition most real devices have a width that is much 
smaller than their length and the junctions are distributed symmetrically so that it 
is reasonable to reduce the problem to one dimension. To do this we expand ip on 
transverse Fourier modes 

(2.11) ¥ ,( a .,„ )t ) = n^_D + ^ n ( M ) C0S (^) , 

n=0 

where j = I/w and the first term takes care of the boundary condition (|2.8[) . After 
inserting (|2.1ip into (|2 .6[) and integrating across y we get the evolution of fay 

(2.12) fa t - cj) xx + —g(x, y = 0)(nfa t + afa + sin0) = v- , 

w I 

where we omitted the in (j)Q and terms in fa , i > 1 which are small because of the 
smailness of the current [17]. 

The boundary conditions are 

(2.13) fa\*=o = H - (1 - v)l, 4> x \ x=l = H + (1 - 



and 



0, elsewhere. 



The factor Wj/w is exactly the "rescaling" of \j(= 1) into A e // = J-^j > 1 due to 

the presence of the lateral passive region [T8] . 

As the area of the junction is reduced the total super-current is reduced and tends 
to zero. Small junctions where the phase variation can be neglected and that have 
a significant supercurrent we introduce the following Dirac distribution model. First 
define the function gu{x), 



9h{x) 



hlj , hlj 

elsewhere, 
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notice that g\(x) — g(x). When the junction widths Wj <C Aj we can further reduce 
the problem by neglecting the variation of the phase inside the junctions. We then 
obtain the Dirac delta function distributed model jT51 US] by making h tend to 0. For 
an n junctions device, 

lim g h (x) = -y^O" - a j)- 
Finally we obtain the 8 — ID model for the device, 

n 

(2.14) 4> tt - <f) xx + djS(x - aj)(K<f>u + sin0 + a<f> t ) = v—, 



where 

(2.15) d 



Wjlj 



and where the boundary conditions are given by eq. (l2.13[) . This is the main model of 
the article and we will analyze it in detail. 

3. Preliminary analysis. From the Josephson equation (jl.ip it can be seen that 
4>t is a voltage. In experiments this instantaneous voltage is of very high frequency 
(~ 500 GHz) and can only be detected by making it beat with a well-known source. 
On the other hand the time average voltage can be measured in a fairly standard way. 
This current voltage relation (the I — V curve) is a characterization of the device. It 
is therefore important for the analysis to explain it. We now give some important 
symmetries of the I — V curve. The first one is the periodicity with respect to the 
magnetic field H. This is similar to the one obtained in the static case [9]. 

3.1. Periodicity of the I — V curve with H. The I — V curve of the device 
modeled by eq. (|2.14p and (|2.13p depends on the magnetic field H . We denote I — V\h 
the I — V curve of the device for the magnetic field H . Let us introduce lj = aj+i — aj 
the distance between two consecutive junctions. Let l m i n be the smallest distance lj. 
We define a harmonic array as a circuit where li is a multiple of l m i n for all i. 

Proposition, Periodicity of the device:. For a harmonic circuit, the / — V\h curve 
is periodic with a period 2-k jl m i n . 

Proof:. Let ^ be a solution of (|2. 14[) for a current 7 and a magnetic field H. We 
introduce f{x) = {2tt /l m in){x — ai) and ip(x;t) — 4>{x]t) + f(x). So ip verifies 

n ^ 

(3.1) <j)tt ~(f>xx + y^ J ~ J [ ls ( x _ a j)( K 0*t + sin^ - /) + a <t>t) = v-, 

J=l 

with ip x (0; t)=H + 2jr/k - (1 - v)j/2, and ip x (l; t) = H + 2ir/l t + (1 - vh/2. Since, 
/(a,) = 2fc7r, Vi S {1; n}, then ip is a solution of (|2.14l) for H + H p = H + 2^jl min 
and for the same 7. 

Conversely, by subtracting / from a solution associated to H + H p and a current 
7, we obtain a solution for H and the same current 7. We have shown that for a given 
7 and H, if there is a solution, we can find another for H + 2k-K jl m i n (k G Z). We 
obtain the same I — V curves. 

(3.2) I-V\ h+h =I-V\h. 
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with H p = 2n/l min . 

In the non harmonic case, if the junctions are set such that lj = Pj/qj, where Pj 
and qj are integers, prime with each other, then / — V\h is periodic with period H p 
such that 

( o q\ „ ^ LCM{q l ;...-q n _ l ) 

(3-3) h p = 27r T777TT7 7, 

HCF{p 1 ;...;p n - 1 ) 

where LCM is the Lowest Common Multiple and HCF the Highest Common Factor. 
To prove this write f(x) = H p (x — ai) and use again the previous argument. 

3.2. The high voltage approximation . When the voltage <f>t is large, the 
phase <p is rotating fast so that one can write 

(3.4) <j>(x,t) = Vt + i/;(x,t), 
where the average 

1 f t+T 

(3.5) ^=TJ ^,t')dt' ^^{x). 
Plugging the ansatz (|3.4p into (|2.14p and taking the average we get 

n 

(iptt) - 4>v xx + d i S ( x - a i) 
i=i 

(3.6) 

[k (iJju) + (sin(W) cos(V')) + (cos(W) sin(V')) + a (ip t ) + aV] = 7. 

Then if we neglect the nonlinear terms we obtain the static equation, we obtain a new 
equation such that {ip t t) — (ipt) — is a solution. Thus, 

n 

(3.7) -(f> Vxx + aV^2d j S(x-a j ) = u-j, 
together with the boundary conditions 

(3.8) cb vx \ x=0 = H - (1 - ^> vx \ x=l =H+{l-v) 1 -. 



By integrating the equation (|3.7[) one sees that this problem has a solution if 

7 



(3.9) V = 



Let us write this high voltage solution for a device with many junctions. At each 
junction x = dj, the phase must be continuous. Integrating equation (|3.7|) over a 
neighborhood of Oj and taking the limit of this neighborhood going to zero one gets 
the jump condition 

(3.10) -[faA+aV?f=0. 
Notice also that outside the junctions, the solution has to be of the form 

<t>v{x) = ~ v Yl x2 + ° lX + ° 2 ' 
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It is then natural to build the solution by steps, using a sort of shooting method. To 
simplify the discussion, let us consider a device with two junctions placed respectively 
at the positions a±, a 2 . We write 



(3.11) 4> v {x) -- 

At each junction aj we have <fij + i(a,j) = <f>j(a,j) and 



4>a{x), 0<x<ax, 
<t>i[x), ai < x < a 2 , 
4>2(x), a 2 <x <l. 



, J+lx (a+)-h x (aj) + aV-f=Q. 



7 2 , 
- V 2l X + 



H-{1 



It is then natural to write 
(3.12) Mx) 

From the relations at each junction, we infer that 

di 

<Pi(x) = (po{x) +1— 

(3.13) 



x + C. 



<f>i(x) = O (x) + 7 

02 (x) = 01 (x) + 7 



di + d 2 

d 2 
d\ + d 2 



(x - ax), 
(x - a 2 ). 



It can be checked that the right boundary condition for x — I can verified by this 
formulation whatever the value of H by choosing 7 thanks to cq. (|3.9|) . The function 
0„ (x) defined by the equations Q3. 12113. 13| is then the solution of the problem (|3.7|3.8|) . 

3.3. Simplification of the S — ID model. The static part 0„ of the high 
voltage solution can be used to simplify the formulation of the problem, in particular 
the boundary conditions. For that we introduce 

(3.14) ^{x,t) = <j>(x,t) -0„(aO, 

so that the sine-Gordon equation (|2.14j) becomes 



(3.15)V*t - V> 



^djSix-aj) 



Kiptt + aipt + sin(V> + 0u) 



Efe=i d k 



= 



with the homogeneous Neumann boundary conditions 

(3.16) 1px\x=0 = 0, 1px\x-l=0. 

Equation p. 151) now contains for each junction a sine term with an argument that is 
shifted by 4> v (aj) — <f> v {a\). When averaging over time we get 

(0( a j) - 0(ai)) t « 0„(%-) - <j> v {a\) ■ 

We define as the phase-shift of the junction j 

(3.17) = <t>v{aj) -^,(01) . 
By a simple translation, equation (|3.15l) becomes 



(3.18)^ tt - ip xx + ^ djS(a 
3=1 



Efe=i d k 



= . 



The boundary conditions (|3.16[) are unchanged and "fi = 0. 
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4. Resonances for an array of point junctions. 

4.1. Spectral problem. When the system is on a resonance, the solution is 
periodic so that it can be written as ip(x,t) = e lut Lp(x). In that case the terms 
aipt + 7/ y\_ 1 dk globally balance each other. The sine term can be averaged out. 
Therefore to satisfy the equation (|3.18[) it is necessary that 



(4.1) 



3=1 







together with the homogeneous Neumann boundary conditions (13. 16[) . 

When the capacities per unit area of the junctions and passive regions are equal 
K = 0, so the array resonates on a cosine Fourier mode, see [T51 119] , 



ip n (x,t) 



' cos 



\~~r) 



In this description, we have neglected the higher harmonics which decay exponentially 
as shown by numerical calculations. When k ^ 0, the eigenmodes of equation (|4.1I) 



differ from nn/l. To analyze them we substitute ip(x,t) 
obtain the following eigenvalue problem 



e iut Lp(x) into (|4~T|t and 



(4.2) 



1 + ^2 d j S{x - a J 



<p = 0, 



together with homogeneous Neumann boundary conditions for (p. 

We consider the two junctions case. We will obtain results that can be generalized 



to an array with n > 2 junctions. We introduce for ease of notation Kj 
then natural to assume a solution 



ndj. It is 



(4.3) 



<p(x) = 



Ax cos{(jjx), 
Ai cos(wx) - 
A 3 cos(wir) - 



B2 sin(wx), 
i?3 sin(wx), 



for 
for 
for 



< x < a\, 



< 
ax L 

d2 < X < I 



ax < x < a,2, 



where the form of <p in the first interval was chosen to satisfy the boundary condition 



at x = 0. As usual <p is continuous at the junctions x 
jump of the derivative is 



a, and one can see that the 



(4.4) 



r 1 2 

[ip x \ i + KjUi <p{a.j) = 0. 



Writing these In — 4 conditions at the junctions and the boundary condition at x = /, 
one obtains a 5th order homogeneous system in Ax, A2, B2, A3, B3. The solution is 
non trivial if the following determinant is zero 



(4.5) 



Cx 


-Cx 


-Si 








Sx 


—Si + kxujCx 


Cx + KxijjSx 











c 2 


s 2 


-C2 


-s 2 





s 2 


-C2 


—S 2 + C2 - 


h K2C0S2 











-St 


Q 



0. 



where Ci = cos(ojai), 5i = sin(a;ai), C'2 — cos(a;a2), ^2 = sin(wa2), Ci — cos(u>l), Si = 
sin(u)l). When more junctions are present in the device, the determinant giving the 
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Fig. 4.1. Dependence of the eigenvalues (zeros of the dispersion relation \4- 7]) ) on the position 
of the junction a\ in the microstrip (left panel) and on the capacity miss-match k (right panel). For 
the left panel, k = 0.5. For the right panel ai = 0.31 = 0.3tt. 



dispersion relation can be generated by adding the elementary component in rows 3 
and 4 corresponding to each additional junction. 

We now give the resonant frequencies when there is only one junction. The 
determinant (|4.5[) becomes 



(4.6) 

which implies 
(4.7) 



C\ —C\ —Si 

Si —Si + KiUjCi Ci + KitvSi 

-Si Ci 



sin(wZ) + nujdi cos(wai) cos(ui(l — ai)) = 0. 



Several remarks should be made on this relation. 

• First, when n <C 1 we recover the usual sin(wZ) = dispersion leading to 
harmonic frequencies. 

• The opposite limit k ^ 1 is more interesting because it leads to a splitting 
of the oscillations in the left and in the right side of the film. We obtain 
cos(ujai) = or cos(uj(1 — ai)) — leading to wai = (2n + l)n/2 or uj(1 — a t ) — 
(2m + l)7r/2 where m, n are integers. 

Notice that Larsen et al [3D] considered a centered junction «i = !/2 in a mi- 
crostrip. Their dispersion relation 

/ Ci f 

tanTr— = -7T— — , 

Je Je 

is similar to the one we get with our approach except that the coefficient is different. 
We obtain using (|4.7p 



(4.8) 



(j. 

tan(u^) = — (jjl ( — f — 1 



3 

wl 



This latter expression gives the correct eigenmodes as shown in the current voltage 
characteristics computed numerically shown in the next sections. 
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4.2. Eigenvectors and inner product. We introduce here the inner product 
associated with the dispersion relation (|4.7p . The eigenvectors ipi (respectively ipj) 
associated to the eigenvalue u>i (resp. u>j) satisfy 

(4.9) <p ixx + ufcpi (1 + S(x - oi)«n) = 0, 

(4.10) <p jxx + ulcpj (1 + 5{x - ai)«i) = 0, 

together with the boundary conditions 

( 4 - u ) fixlxeW] = °' and ^*L{o,i} = °- 

We assume the eigenvalues to be different. As usual we multiply eq. (|4.9[) by ifj and 
eq. (|4.10p by <pi , substract the second from the first and integrate on the domain. We 
obtain 



(4.12) / ((p ixx <Pj - tp Jxx ipi)dx + (w, - ujj) 



i 

(PiPjdx + Kiipi(ai)ipj(ai) 



= 0. 



The first integral is zero because of the boundary conditions (|4. 1 If) . Since u>i ^ u)j we 
get 

(4.13) / (1 + k\S(x — ai))(fiifjdx = 0. 

Jo 

Thus we obtain the orthogonality of the eigenvectors (fi, i is an integer with the 
associated inner product defined by 

(4.14) (f;g}= [ (1 + k x 5{x - ax)) fgdx. 

Jo 

It is then possible to normalize the eigenvectors so that they form an orthonormal 
basis. When there are n junctions in the array, the inner product can be generalized 
easily to 



(4-15) (f-g) 



J ^1 + ^ K t 8{x - aj^j fgdz 



In the case of a single junction, the normalized eigenvectors (p n (x) are given by 

J A n cos(uj n x), for 0<x<ax, 

{ } ^ X >-{ A n J^ a \ )) cos(u Jn (l-x)), for ai <x<l, 

where to n satisfies the dispersion relation (|4.7I) and 
(4.17) A n = ' 



I | sii^o^ai i cos 2 u„ai sin 2u„ (Z-ai ) . 2 , , „ 

77 ~r -* i 3 77 7 ~ A ~r ™1 COS UJ 71 Cll 

Fig. 14.21 shows the 4 non trivial eigenmodes 4>i, i = 1 — 4 for a large capacity miss- 
match. Notice how the modes are almost zero on one side of the cavity. 
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0.3n 7i 0.3 jc 7i 

x x 

Fig. 4.2. Normalized eigenmodes ipi(x), i = 1, 2 (top panels from left to right) and <j>i( x )> « = 3, 4 
(top panels from left to right) for a large capacity miss-match k = 30. The other parameters are 
dx = 1, I = 7r and a\ = 0.3-7T. 

4.3. Eigenvalues and I — V curves. The eigenvalues calculated in the previous 
section appear as resonances in the current-voltage characteristics of the device. We 
computed the I—V curves from the numerical solution of equation (|2.14[) . The singular 
partial differential equation is integrated over reference volumes and the time advance 
is done by an ordinary differential equation solver (see Appendix [5] for more details) . 
The average voltage V = ((fit) is computed over a time interval 10/a after waiting a 
time 100/a for the solution to stabilize. When the system is locked on a resonance, it 
oscillates periodically on an eigenfrequency solution of the dispersion relation eq. (14.71) . 
We will show that the numerical solution follows closely the corresponding eigenvector. 
We choose a device with I = ir, a± = 1/3, d\ = 1 and take three limiting cases n = 
K = 1 and k = 30. When there is no capacity miss-match (n — 0) we showed in [TS] 
that the resonances of the I — V curve are positioned at V — kn/l with k integer and 
arc bounded by 

7? = diaV, and 7/ = dx^/iaVY^T. 

This is shown in the top panel of Fig. 14.31 In the middle panel we show the I — V 
curve for k = 1. Notice how the resonances are not equally spaced and get sharper. 
The resonances are very close to the ones given by the dispersion relation, they are 
reported in table (|4.3[) . Finally on the bottom panel we computed the I — V curve 
for k = 30. This large value is typical in experiments where the oxide layer in the 
junction is about 100 times thinner than in the passive region. This together with 
the ratio d\ = Wj/w gives k\ ~ 50. Notice how the resonances are vertical indicating 
that the system is almost linear. This is to be expected because k » 1 so that 
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i 




cos(wia) 


Normalization Ai 





0. 


1. 


0.491 


1 


0.917 


0.649 


0.562 


2 


1.931 


-0.247 


1.094 


3 


2.438 


-0.665 


0.704 


4 


3.682 


-0.946 


0.237 



Table 4.1 

Eigenvalues for the array of Fig. \4-3\ with k = 1 (middle panel of Fig. \4-3\ l- 



the sine term plays little role. The position of the resonances is exactly given by 
the dispersion relation (|4.7p . They are reported in the table (|4.3p together with the 
coupling coefficient <p n (ai). The eigenvalues are clearly not harmonic. Notice how 
some resonances do not go all the way to the junction curve. Another observation is 
that we have a resonance for w 2 even though the coupling coefficient is almost zero. 
We will come back to this point in the next section. Also it is interesting that there 
is no resonance for uj^. 

Fig. 14.41 shows the eigenvectors in dashed line (green online) for a realistic value 
K\ = 30 corresponding to a capacity miss-match k — 30 and a surface ratio d\ = 
w\ h/w = 1. On the same plot are indicated in continuous line (red online) the 
instantaneous voltages tpt for different successive times. From top to bottom panels 
we show respectively the first, second and third resonance. One can see that the 
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i 




cos(ojia) 


Normalization Ai 





0. 


1. 


0.1737 


1 


0.734 


0.770 


5.341 IQ-' 2 


2 


1.687 


-1.950 IQ-' A 


1.447 


3 


2.150 


-0.440 


3.462 10~ 2 


4 


3.575 


-0.974 


9.102 10~ a 



Table 4.2 

Eigenvalues for the array of Fig. \j.3\ with ft = 30 (bottom panel of Fig. 




voltages follows very well the eigenvectors. 

4.4. Projection on the normal modes. The normal modes that we have ex- 
hibited can be used to project the solution. This enables to analyze the dynamical 
behavior in a simple way. The eigenvectors solutions of the problem (|4.9I) are or- 
thonormal for the inner product defined above. We will assume that all solutions </> 
belong to the vector space generated by the eigenvectors ifii. This is difficult to justify 
theoretically but we will come back to this using the numerical results in the last 
section. Then the solution ip of equation (|3.15p <f> can be written as 



(4.18) 



+ 00 

i=0 



14 



J.G. Caputo and L. Loukitch 



We assume the uniform convergence of the series (|4. 18[) so that we can permute 
integrals and sums. 

We will present the calculations for a single junction for simplicity. The results 
can be generalized to arrays with multiple junctions. We replace ip in eq. (|3.15p for 
one junction, 

+00 

(4.19) 1=1 

di5(x - d) f i^PiVi + uP'iVi) + sin (^ + #1) - j- J = 0, 

We multiply ea. (|4.19l) by fj{x) and integrate it on its domain 

el \ fl 



+00 / .1 

<=i V Jo 



^^•cte + Ki^dVjCd) J / (p"ipjdx + dia/3' i ip i (ai)cpj(ai) 



=0 if i^J, =1 else 

(4.20) +disin(V>(ai) + *i)^(ai) -7^(01) = 0. 
and from eq. (l4.2[) we know that, 

If" = - (1 + K\5{X - Ol)) W, 2 (y5i. 

Eq. (|4.20p becomes 

+00 , 

/3" + ^/3jW? / + _ ai))<pnpjdx-\ — — /3 t Vi(dVj(d) 



i=l 



-0 if —1 else 



(4.21) +y sin (^(oi) + *i)^-(oi)- 7^(01) =0. 
We obtain the final equation giving the evolution of /3j in terms of <f> 

(4.22) + wfa + Cj (a<f> t (ai) + sin (0(d)) - = 0, 

where we have returned to the usual field <f> and where the coupling coefficient Cj is 

(4.23) Cj = diipj(a\). 

This coefficient does not depend directly on k. Also notice that all equations are 
coupled by the same term 

(4.24) F = a<t> t (ai) + sin(0(d)) ~ ]~, 

di 

where the coefficient Cj regulates the forcing for each mode. When n junctions are 
present in the device, the modal equations can be generalized to 



n , 

(4.25) % + u% + E C M a<t)t ( flfc ) + sin ^h)) - Y^T" 

k=i v ^=1 ai 

where the generalized coupling coefficient is 

(4.26) d=d k <pj{a k ). 



= 0, 
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Fig. 5.1. Plot of the amplitudes 0i(t), i = 1 — 4 obtained by projecting the numerical solution 
for k = 1 at the top of the second resonance onto the normal mode <pj . The modes are fi\ continuous 
line, (red online), 02 dashed line (green online), 03 short dash (blue online) and 04 dotted line (pink 
online). In the left panel the projection is done directly on the solution <j>(x,t). In the right panel 
we have subtracted the high voltage solution <j> v (x) for the given current 7 = 1.36. The time unit is 
the period T2 = 2ir/uJ2 ~ 3.25. 



5. Numerical analysis of the IV curves. To analyze the mechanism leading 
to a resonance in the IV curves shown in Fig. 14.31 we project the numerical solution 
onto the normal modes that we defined in section 4. Following the definition of the 
inner product (|4. 14[) we have 

ft = ((j>\(j)i) = / <p4>idx + Kdi(f>(ai)<f>i(ai). 
Jo 

The integral on the right hand side is computed using the trapeze method. Fig. 15.11 
shows a plot of the amplitudes ft, i = 1 — 4 for the second resonance with k = 1 
(middle panel for Fig. I4.3[) . Clearly ft is dominant. In the left panel we did not 
subtract the high voltage solution <j> v (x) Q3. 12113.13]) so that the other modes appear 
as parasites. If the high voltage solution is taken out then we have a clear dominance 
of ft, all the other modes being close to 0. This shows that we have to a good 
approximation 

(5.1) <f>{x, t) = i/j(x, t) + 4> v (x) Po(t)<p (x) + Pi(t)<Pi(x) + <t> v {x), 

where i = 2 and where we have included the mode that is always present. We 
recover the results suggested by the plots of Fig. 14.41 We have observed this for all 
the resonances in the I — V curve. For example for k = 30 (bottom panel of Fig. I4.3[) 
we show in Fig. 15.21 the amplitudes ft from top to bottom for the 1st , 2nd and 3rd 
resonance. Again the dominant amplitudes are from top to bottom ft, ft and ft. 

Therefore the solution at the top of the resonances is given by (|5.1j) to a good 
approximation. To analyze how we reach this state we have projected the solution 
on the normal modes for increasing values of the current 7 all the way to the top 
of the resonance. The calculations were done over a long time interval (about 200 
periods) . Projecting the solution tj> we observed a drift in the amplitudes ft due to the 
rapid increase of <fi and the finite precision of the evaluation of the integrals. To avoid 
these technical problems we have projected the time derivative cj>f The qualitative 
conclusions are the same as for except that we will look at ft t . Fig. 15.21 shows 
three amplitudes as a function of time for 7 = 1.12, 1.3, 1.42 and 1.56 and a voltage 
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Fig. 5.2. Plot of the amplitudes ft(t), i = 1 — 3 obtained by projecting the numerical solution 
at the top of the second resonance onto the normal mode ifii for re = 30. The dominant mode is 
resp. i = 1,2,3 for resp. the 1st resonance (top panel), the second resonance (middle panel) and 
the third resonance (bottom panel). 



V ~ L03 near the 3rd resonance for n = 1. Only three periods T3 = 2tt/uj^, have been 
represented for clarity, the rest of the time evolution is the same. In the top left 
panel for 7 = 1.12, the amplitude of the mode 3 is about 1 with small components 
in the modes 2 and 1. When the current is increased the amplitude of the 3rd mode 
increases and becomes periodic of period T3. The other modes tend rapidly to 0. 

It is instructive to compute numerically the forcing term F as one progresses 
up the resonance. Fig. 15.41 shows F(t) for three periods T3 for the four values of 7 
analyzed in Fig. 15.31 The amplitude of F decreases for increasing 7 and F becomes 
periodic of period T3. This explains why we obtain the correct resonant modes using 
the spectral problem (|4.2[) . Note that when n = the forcing term F of the amplitude 
equations tends to when one gets close to the top of the resonances. Then one can 
solve the differential equation for <f>{a\ \ t) 

(5.2) acj>' '(ai; t) + sin ( 0(cn; t)) -^-=0, 

and close the system by expanding this solution using the standard Fourier modes [15]. 
This is not the case when k ^ as shown by these numerical results. To analyze the 
resonances, we assume as in [12] that when in resonance, the solution has the spatial 
structure of the corresponding cigenmode. 

(5.3) ip(x,t) = /3 (t)cpa + /3 n (t)tp n (x), 



where the first term on the right corresponds to the zero mode. The evolution of 
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Fig. 5.3. Plot of the amplitudes of 4>t(x,t) as a function of time for four different values of 
the current, from top left to bottom right 7 = 1.12, 1.3, 1.42, 1.56, on the third resonance V = 0J3 for 
K = 1. The index of the modes are i = 1 continuous (red online), i = 2 dashed line (green online), 
i = 3 short dash (blue online) and i = 4 dotted line (pink online). The time interval is the period 
T 3 = 2tt/uj 3 ss 2.58. 
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Fig. 5.4. Plot of the forcing term F from \4-2J$ as a function of time for four different values 
of the current on the third resonance V = L03 for k = 1. The values are 7 = 1.12 continuous (red 
online), 7 = 1.3 dashed line (green online), 7 = 1.42 short dash (blue online) and 7 = 1.56 dotted 
line (pink online). The time interval is the period T3 = 2tt/lJ3 pb 2.58. 



Po, p n is then given by equation (|4.22l) 



(5.4) 
(5.5) 



PI + oj 2 n P n + c n F = 0, 
$ + c F = 0, 
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Fig. 5.5. Plot of the frequency ui vs. m for a device with two junctions where ai = 0.37T, a 2 = 
0.5n, di = d 2 = 0.1. 

where the forcing term is 

(5.6) F = a/3 n ip n + a/3' ip + sin (p n (p n + /3 Q (p ) - — 

«i 

To understand the specific shape of the resonances, the fact that we cannot obtain in 
the IV curve the right part of the resonance curve, one could carry out a bifurcation 
analysis similar to the one of [21]. However this is out of the scope of this article. 

The situation is more complex when there are more junctions in the device. As 
an example we consider a two junction device with I = tt, a\ — 0.37T, a-i = 0.57T, d\ — 
d 2 =0.1. The eigenmodes oj n are plotted in Fig. I5.5l as a function of k and one can see 
them shift from integer values even for small k. In Fig. 15.61 we plot the I — V curves 
obtained for this device for k — 0, 4 and 8 with d\ =0.1 so that n\ = d\K = 0, 0.4 
and 0.8. The resonances observed correspond to the eigenfrequencies obtained. For 
k = we have explained the height of the resonances using an approximate theory [19] 
based on the amplitude of oscillation of 4>t(o-j) for each junction. This is the envelope 
function plotted in dashed line in the top left panel of Fig. 15.61 When n > it is 
more difficult. We think that this amplitude of oscillation, which can be obtain by 
the eigenvectors does not determine completely the height of the resonances. When 
the length of the device is larger more resonances can be accommodated in the I — V 
curve. An interesting effect we found is that for large values of k the system locks to 
linear combinations of the eigenfrequencies. Such an example is shown in Fig. [S] for 
a two junction device in a microstrip of length / = 10. For k > 5 there appears in 
the I — V curve a resonance for u>i + ujq. For k = 10 (bottom left panel) we see in 
addition resonances for 2lj\, uj\ + ujq. This is typical of a linear system. When we 
observe numerically the evolution of 4>t at this resonance we see the two eigenvectors 
1 and 9. But we cannot explain this linear behavior: how the system can sum two 
non linear solution? 

5.1. Conclusion. To summarize we have analyzed a long wave model describing 
a parallel array of Josephson junctions. We defined an appropriate spectral problem 
whose spectrum gives the resonances of any array with junctions of arbitrary (small) 
sizes and positions. This task would not have been possible without the complexity 
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Fig. 5.6. IV curves for the two junction circuit studied in Fig. 15.51 The parameters are 
I = 7r, d\ = d2 = 0.1, a = 0.25, ft = 0. From top left to bottom k = 0, 4 and S. Jra t/ie top panel the 
envelope function (see text) is plotted as a dashed line (pink online). 



reduction provided by asymptotic analysis. 

The adapted spectral problem leads to an inner product so that it becomes pos- 
sible to project the dynamics of the system and describe arrays with more junctions. 

It may now be possible to solve the inverse problem of finding the device yielding 
a given I-V curve. Another open question is the study of the amplitude equations 
(|5.4|) to analyze the stability of the resonances. 
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6. Appendix. The basis of the method is to discretize the spatial part of the 
operator and keep the temporal part as such. We thereby transform the partial differ- 
ential equation into a system of ordinary differential equations. This method allows 
to increase the precision of the approximation in time and space independently and 
easily. In our case the operator is a distribution so that the natural way to give it 
meaning is to integrate it over a volume. We therefore choose as space discretisation 
the finite volume approximation where the operator is integrated over reference vol- 
umes. The value of the function is assumed constant in each volume. As solver for the 
system of differential equations, we use the Runge-Kutta method of order 4-5 intro- 
duced by Dormand and Prince implemented as the Fortran code DOPRI5 by Hairer 
and Norsett [22] which enables to control the local error by varying the time-step. 

We first transform (..) into a system of first order partial differential equations 
We write ip(x,t) = dj t (x,t). 



S(x — a)(Kipt{x, t) + cnp(x, t)) + di sm(4>(x, t))) + vj/l 



ip(x,t) = <pt(x,t) 
ipt{x,t) = 4> xx {x,t) 
(6.1) 

with the boundary conditions : cj) x \i = H — (1 — is)j/2, and (f> x \_i = H+ (1 — u)^/2. 

For simplicity we will describe the implementation of the finite volume discretisa- 
tion in the case of a single junction. We introduce reference volumes Vk whose centers 
we call Xk, 1 < k < nn. The discretisation points are placed such that the point 
Xng+i is at the junction, (x ng +i = a). We thus define Xk and Vk using the following 
identities 

h a 



Vk = 



Xk 



Kg 

— ,Xk 

2 



with (ng + l)h g = a 



h d 

x k - -^,x k 



hd 
2 



< k < ng + 1 



, ng + 1 < k < nn + 1 



with (nn — ng)hd — I — a. Finally at the junction, k = ng + 1 



V K 



j •Eng-\-l 



hd 
2 



nn, ng and nd are respectively the total number of discretisation points, the number 
of points to the left of the junction and the number of points to the right. 

For a fixed t, we assume <j>(x, t) to be constant on each volume Vfc, so that 



Xh + - 



(x, t)dx = h(f>(xk,t) , with h = hg or h = hd 



Integrating over Vk yields: 

In the linear part of the partial differential equation : < k < nn + 1 and 
k 7^ ng + 1 : 



(6.2) 



ip(x k ,t) = (j> t (x k ,t) 
ip t (xk t) = ( t>( x i'+-i-> t )-' 2 <t>( x k,t)+4>(xk-i,t) 



+ .] 
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with h = hg for < k < ng + 1 or h = hd for k > ng + 1 . We recognize the usual 
discretisation of the second derivative. 

At the junction: k = ng + 1, we obtain 



I 



5(x - a) (n(j)tt(x,t) + a4>t(x,t) + d\ sm(4>(x , t))) 



dt sm(<f>(x ng+1 ,t)) + a4> t (x ng+1 ,t) 

So that the final system is: 
ip{x ng+ i,t) = 4> t (x ng+ i,t) 

, , ,n _ _4 / <t>{x ng +2,t) ~ <j>(x ng+ i,t) _ (j)(x ng+1 ,t) - 4>{x ng ,t) 

mXn3+1,) ~[hg + hd\ hg/2 hd/2 

2 1 1 

(di sm(<t>(x ng+1 ,t)) + a(j} t (x ng+1 ,t)) +j 



hg + hd 



1 — rffeK 



The previous system of ordinary differential equations is then integrated numeri- 
cally using the DOPRI5 integrator. 



